Métodos Clássicos de Descida com Busca de Armijo¶

Otimização Não Linear¶

Funções Auxiliares¶

In [1]:
# importando funçõs necessárias
!pip install numdifftools
from IPython.display import clear_output
import numpy as np
import pandas as pd 
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from matplotlib import cm
import numdifftools as nd
import warnings
clear_output()
In [2]:
pd.set_option('display.float_format', lambda x: '%.8f' % x)
np.set_printoptions(8)
warnings.filterwarnings("ignore")
In [3]:
# funções auxiliares
def plot_function(f, title, ndim, dom = np.linspace(-5, 5, 500), angle = (20, 20)):

  """
    Função que plota gráficos em 2-d ou 3-d.

    Parâmetros
    ----------------------------------------
    f : callable
        Função objetivo (função custo).
    title : str
        Título do gráfico
    ndim : str (2 ou 3)
        Dimensão da função
    dom : list or np.arange or np.linspace, opcional
        Domínio da função desejado 
    angle : tuple,  opcional
        Ângulo para ver o gráfico. Padrão = (20, 20)
    Saída
    ----------------------------------------
    plot
  """
  
  
  plt.style.use('fivethirtyeight')
  
  # 2-d plot: y = f(x)
  if ndim == 2:

    X = dom     
    Y  = f(X)

    fig = plt.figure(figsize=(8,8))
    ax = fig.gca(projection='2d')
    ax.set_title(title)
    ax.set_xlabel('$x$')
    ax.set_ylabel('$f(x)$')
    ax.plot(X, Y)
    ax.view_init(angle[0], angle[1])
    plt.tight_layout()
    plt.show()

    
    #3-d plot: z = f(x, y)
  else:

    x, y = dom, dom 
    X, Y = np.meshgrid(x, y)
    Z = f([X, Y])
    fig = plt.figure(figsize=(8,8))
    ax = fig.gca(projection='3d')
    ax.set_title(title)
    ax.set_xlabel('$x_1$')
    ax.set_ylabel('$x_2$')
    ax.set_zlabel('$f(x_1, x_2)$')
    ax.plot_surface(X, Y, Z, cmap='jet')
    ax.view_init(angle[0], angle[1])
    plt.tight_layout()
    plt.show()


def plot_results(f, xi, yi, method, dim = [np.linspace(-7, 7, 500), np.linspace(-7, 7, 500)]):

    """
    Função que plota os result

    Parâmetros
    ----------------------------------------
    f : callable
        Função objetivo (função custo).
    xi : np.array
        Valores de x em cada iteração do método
    yi : np.array
        Valores de y em cada iteração do método
    method: str
        Método desejado: '1', '2' ou '3'
    dim : np.linspace, opcional
        Valores de x e y para para plotar as curvas de nível 
    Saída
    ----------------------------------------
    plot
  """


    plt.style.use('fivethirtyeight')
  
    if method.lower().strip() == '1': title = 'Descida do gradiente com busca de Armijo' 
    elif method.lower().strip() == '2': title = 'Método de Newton com busca de Armijo'
    else: title = 'Método BFGS com busca de Armijo'

    # z = f(x,y)
    if xi.shape[1] == 2:
      
      x, y = dim[0], dim[1]
      X, Y = np.meshgrid(x, y)
      Z = f([X, Y])
      fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 8))
      plt.suptitle(title, y = 1.05)

      # curvas de nível
      ax1.plot(xi[:,0], xi[:,1], linestyle='--', marker='o', color='black', linewidth = 3)
      ax1.plot(xi[-1,0], xi[-1,1], 'ro', markersize = 11)
      ax1.set(title='Caminho durante a otimização - Curvas de Nível', xlabel='x1', ylabel='x2')
      CS = ax1.contour(X, Y, Z, 15, cmap = 'jet')
      ax1.clabel(CS, fontsize='smaller', fmt='%1.2f')
      
      # valor da função custo em cada iteração
      ax2.plot(yi, linestyle='--', marker='o', color='black')
      ax2.plot(len(yi)-1, yi[-1], 'ro', markersize = 11)
      ax2.xaxis.set_major_locator(MaxNLocator(integer=True))
      ax2.set(title = 'Valor da função custo durante a otimização', xlabel='Iterações', ylabel='Valor da função custo')
      ax2.legend(['Busca de Armijo'])
    
      plt.show()
    
    # y = f(x)
    else:

      x = dim[0]
      y = f(x)
      fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 8))
      plt.suptitle(title, y = 1.05)

      ax1.plot(xi[:,0], xi[:,1], linestyle='--', marker='o', color='black', linewidth = 3)
      ax1.plot(xi[-1,0], xi[-1,1], 'ro', markersize = 11)
      ax1.set(title='Caminho durante a otimização', xlabel='x1', ylabel='x2')
    
      ax2.plot(yi, linestyle='--', marker='o', color='black')
      ax2.plot(len(yi)-1, yi[-1], 'ro', markersize = 11)
      ax2.xaxis.set_major_locator(MaxNLocator(integer=True))
      ax2.set(title = 'Valor da função custo durante a otimização', xlabel='Iterações', ylabel='Valor da função custo')
      ax2.legend(['Busca de Armijo'])
    
      plt.tight_layout()
      plt.show()



def printa_resultados(chute, xi, yi, grad_norm):


    chute = np.round(chute, 8)

    print('\n*************** Resultados ***************')
    print()
    print('Número de iterações: ', len(yi))  
    print()
    print('***** Primeira iteração *****')
    print()
    print(f'f({chute}) = {np.round(yi[0], 8)}')
    print()
    print(f"||∇f({chute})|| = ", np.round(grad_norm[0], 8))
    print()
    print()
    print(f'***** Última iteração: {len(yi)}° *****')
    print()
    print(f'f({np.round(xi[-1], 8)}) = {np.round(yi[-1], 8)}')
    print()
    print(f"||∇f({np.round(xi[-1], 8)})|| = ", np.round(grad_norm[-1], 8))


def minimiza_f(f, metodo, chute, sigma = 0.02, tol = 1e-8):

  """ Função para testar os métodos

  Parâmetros
  ----------------------------------------
    f : callable
      função objetivo no formato estabelecido
    
    metodo : int (1, 2 ou 3)
      1 - Método do Gradiente
      2 - Método de Newton
      3 - Método BFGS
    
    chute : np.array ou list
      array com o ponto inicial (mesma dimensão da função)
    
    sigma : float > 0, opcional
      cte de decréscimo de Armijo. Padrão - 0.02
    
    tol : float > 0, opcional
      tolerância pré-definida. Padrão = 1e-8 

  Saída
  ----------------------------------------
    Resultados obtidos da minimização """


  f_grad = nd.Gradient(f)

  m = {1:'Método de Descida do Gradiente com condição de Armijo', 2:'Método de Newton com condição de Armijo', 3: 'Método Broyden–Fletcher–Goldfarb–Shanno (BFGS) com condição de Armijo'}
  print(f"{m[metodo]} \n\n")

  if metodo == 1:
    xi, yi, grad_norm =  GradientDescent(f, f_grad, chute = chute, sigma = sigma, tol=tol)
    printa_resultados(chute, xi, yi, grad_norm)

  elif metodo == 2:
    f_hessian = nd.Hessian(f)
    xi, yi, grad_norm =  Newton(f, f_grad, f_hessian, chute = chute, sigma = sigma, tol=tol)

    if len(xi) != 0:
      printa_resultados(chute, xi, yi, grad_norm)

  else:
    xi, yi, grad_norm =  BFGS(f, f_grad, chute = chute, sigma = sigma, tol=tol)
    printa_resultados(chute, xi, yi, grad_norm)
  
  if len(chute) == 2 and len(xi) != 0:
    print('\n\n\n\n')
    plot_results(f, xi, yi, dim = [np.linspace(int(np.min(xi[:,0]))-7, int(np.max(xi[:,0]))+7, 1000), np.linspace(int(np.min(xi[:,1]))-7, int(np.max(xi[:,1]))+7, 1000)], method = str(metodo))
In [ ]:
 

Introdução¶

O Problema da Minimização¶


Uma das principais aplicações do Cálculo Diferencial e Integral é na área de otimização, uma vez que, na maioria das vezes, o objetivo é maximizar ou minimizar determinada função.

Alguns exemplos disso vão desde a otimização de carteira, no mercado financeiro; caminhos, para serviços de entrega e sistemas GPS; e em diversos modelos de aprendizado de máquina, nos quais a função custo deve ser minimizada para se obter os melhores parâmetros que representarão adequadamente seus dados para determinada tarefa.

Nisso, a otimização não linear apresenta várias maneiras de lidar com tais problemas. Neste trabalho serão apresentados os métodos clássicos de descida, sendo eles: Método do gradiente, de Newton e BFGS, todos com busca linear de Armijo. Junto com cada método será apresentada sua descrição teórica, de algoritmo e uma breve análise da minimização de algumas funções, que foram retiradas de exemplos dados em $[1]$, $[2]$ e $[3]$. A parte teórica foi inspirada em notas de aula, em $[4]$ e $[5]$.

A linguagem de programação escolhida foi Python, versão 3.7.13.

In [4]:
!python --version
Python 3.9.7

Algoritmos¶

Algoritmo de Busca Linear com Condição de Armijo¶

A Condição de Armijo (1° Condição de Wolfe) é definida da seguinte maneira:

$$ $$$$ \begin{equation}\tag{1} f(x_{k} + \lambda_{k}d_{k}) < f(x_{k}) + \sigma\lambda_{k}\nabla^{T}{f(x_k)}d_{k} \end{equation} $$$$ $$

Onde:

$x_k \equiv \text{chute inicial (initial guess)}$

$\sigma \in (0, 1) \equiv \text{constante de decréscimo}$

$\lambda_{k} > 0 \equiv \text{passo na iteração } k$

$\nabla^{T}{f(x_{k})}\,d_{k} < 0 \equiv \text{direção de descida em } x_{k}$

$$ $$

Assim, a condição estabelece uma maneira de calcular o passo necessário para um Algoritmo de Descida, como os que serão apresentados abaixo. Além disso, vale ressaltar que este decréscimo da função é proporcional ao passo. Dessa maneira, esquematizando o Algoritmo de busca com condição de Armijo, temos:

$$ $$



Algoritmo 1 - Busca Linear com Condição de Armijo




Dada uma direção de descida $d_{k}$:

$$ $$
  1. Começar com $\gamma = 0.5$ fixo e $\lambda = 1$;
  1. Se a Equação (1) for satisfeita: Fazer $\lambda_k = \lambda$ e $x^{k+1} = x^{k} + \lambda_{k}d_{k}$;
  1. Caso contrário: Atualizar o valor de ${\lambda}$. Fazer $\lambda = \gamma \cdot \lambda$ e ir para o Passo 2;




Importante ressaltar que $\gamma$ é uma constante em que $\gamma \in (0, 1)$.

Abaixo está a implementação do Algoritmo de busca com condição de Armijo (Algoritmo 1).

Algoritmo¶

In [5]:
def Armijo_Search(f, xk, dk, sigma = 0.02, gamma = 0.5):
  
  '''
    Busca linear com Condição de Armijo (1° Condição de Wolfe) - Desigualdade 1, que consiste em fazer um decréscimo da f proporcional ao tamanho do passo 
    
    Parâmetros
    ----------------------------------------
    f : callable
        Função objetivo (função custo).
    xk : array
        Ponto atual.
    dk : array
        Direção de descida (grad_T(f(x_{k})) . d_{k}  < 0).
    sigma : float, opcional
        Valor de sigma entre (0, 1) - constante de decréscimo. Padrão = 0.02
    gamma : float, opcional
        Valor de gamma entre (0, 1). Padrão = 0.5

    Saída
    ----------------------------------------
    lambda : float
        Valor de lambda que satisfaz a condição de Armijo.
    f_x0 : float
        Valor de f no ponto x_{k+1}.
  '''

  # Começa constante lambda = 1

  lambda_ = 1

  # Direção de descida d = grad_{T}(f(x_{k})) . dk < 0  
  desc_direction = np.dot(nd.Gradient(f)(xk), dk)

  f_x0 = f(xk + lambda_ * dk)

  # Calcula a função em um ponto menor que x_k com passo inicial lambda = 1
  # Condição de Armijo para determinar o tamanho do passo, diminuindo a função custo 
  while (f_x0 > f(xk) + sigma * lambda_ * desc_direction):
    
    lambda_ = gamma * lambda_
    f_x0 = f(xk + lambda_ * dk)
  
  # Retorna o passo lambda e o valor da função
  return lambda_, f_x0 

Método do Gradiente com Busca de Armijo¶

O Método do Gradiente pode ser de dois tipos: Subida ou Descida. O interesse nesse trabalho é o de minimização de funções, logo, o de Descida do Gradiente, uma vez que o vetor gradiente fornece a direção de maior crescimento da função.

Esse método consiste em se obter uma direção $d_k$ tal que $d_k = - \nabla f(\vec{x})$ e a partir disso calcular o novo ponto $x_{k+1}$ em que a norma do gradiente $\lVert \nabla f(x_{k+1}) \rVert$ seja menor e o custo da função $f(x_{k+1})$ diminua. O algoritmo esquematizado está abaixo:



Algoritmo 2 - Método do Gradiente com condição de Armijo




Dado um ponto $x_{0} \in \Re^{n}$ ($\textit{initial guess}$) e para $\nabla f(x_{k}) \neq 0$, realizar os seguintes passos:

1. Calcular a direção de descida $d_{k} = - \nabla f(x_{0})$;

2. Realizar Busca de Armijo e encontrar um $\lambda_k$ que satisfaça a desigualdade 1;

3. Calcular o novo ponto $x_{k+1} = x_{k} + \lambda_{k} d_{k} = x_{k} - \lambda_{k} \nabla f(x_{k})$;




Em que tais etapas são realizadas até algum critério de parada pré-estabelecido. No presente trabalho foram 2 critérios: ou a parada se dava por um alguma tolerância referente à norma do gradiente da função ($\rVert\nabla f(x_{k}) \lVert \leq 10^{-18}$), por exemplo, ou por um número máximo de iterações (500 por $\textit{default}$). Ambos os parâmetros são alteráveis nas funções apresentadas.

A convergência do Método de Descida do Gradiente é linear, ou seja, a sequência {$x_k$} converge para $x^{*}$ quando:

\begin{equation} \lVert x_{k+1} - x^{*} \rVert \leq \epsilon \, \lVert x_{k} - x^{*} \rVert \end{equation}

$\forall \, k \in \mathbb{N}$ e com $\epsilon > 0$.

A chamada para a Função que implementa o Método do Gradiente com Condição de Armijo é:

GradientDescent(f, f_grad, chute, sigma, tol)

Onde:

$$ $$
  • f é a função desejada (tipo: callable)

  • f_grad é o gradiente da função f (tipo: callable)

  • chute é o $x_{0}$, ponto inicial (tipo: np.array ou list)

  • sigma é a constante de decréscimo de Armijo (tipo: float $\in (0, 1)$. Padrão $= 0.02$)

  • tol é a tolerância desejada ($\lVert \nabla f(x_{k+1}) \rVert \leq 10^{-8}$) por exemplo (tipo: float $> 0$. Padrão = $10^{-8}$)

$$ $$

E pode retornar até 3 argumentos que são tipo np.array utilizados para fazer gráficos e outras análises visuais.

$$ $$
  • x são pontos $x_k$ calculados pelo Método

  • y são valores de $f(x_{k})$ calculados pelo Método

  • grad_f são valores de $\lVert \nabla f(x_k) \rVert$ calculados pelo Método

$$ $$

Para suprimir algum retorno, por exemplo, o grad_f, basta fazer:

$$ $$
x, y, _ = GradientDescent(f, f_grad, chute, sigma, tol)
$$ $$

Exemplos de como utilizar tais funções serão fornecidos mais adiante com funções pré-estabelecidas. Abaixo temos a implementação do Algoritmo 2:

Algoritmo¶

In [6]:
def GradientDescent(f, f_grad, chute, sigma = 0.02, tol=1e-8):

  """Algoritmo de Descida do Gradiente com busca linear com condição de Armijo. Esquema:

      I - Direção de descida: d_k := −∇f(x).
      II - Determinação do passo com busca de Armijo.
      III - Obtém o próximo candidato.

      Parâmetros
      --------------------
      f : callable
          Função custo
      f_grad : callable
          Gradiente da função f
      chute : array
          Valor inicial de x ("chute")
      sigma : float, opcional
          Constante de Decréscimo de Armijo. Padrão = 0.02
      tol : float, opcional
          Tolerância padrão de 1e-8

      Saída
      --------------------
      xk : array
          valores de xk
      yk : array
          valores da função custo em cada xk
      grad_f : array
          valores da norma do gradiente da f em xk """
    
  # Valores iniciais de xk, fk e grad_fk
  xk = chute    
  fk = f(xk)
  grad_fk = f_grad(xk)
  grad_fk_norm = np.linalg.norm(grad_fk)
  max_iter = 500

  # Inicializa o número de iterações e a lista para fazer os plots dos valores de x e y
  num_iter = 0
  x_pontos = [xk]
  y_pontos = [fk]
  grad_f = [grad_fk_norm]
  print(f'Chute inicial: y = {fk}, x = {xk} \n')

  # Calcula nova iteração com busca de Armijo
  while (grad_fk_norm > tol and num_iter < max_iter):
      
    # Determina a direção
    pk = -grad_fk

    if (np.linalg.norm(f(xk)) > 1e16) or (grad_fk_norm > 1e16):
      print("\nErro: Overflow\n")
      break 

    # Faz a busca de Armijo, obtendo o passo lambda e a função custo naquele passo
    lambda_, fk = Armijo_Search(f, xk, pk, sigma)
    
    # calcula x_{k+1}
    xk = xk + lambda_ * pk
    grad_fk = f_grad(xk)
    grad_fk_norm = np.linalg.norm(grad_fk)

    # Itera mais uma vez 
    num_iter += 1
    x_pontos.append(xk)
    y_pontos.append(fk)
    grad_f.append(grad_fk_norm)

  # print results
  if num_iter == max_iter:
    print('\nNúmero de iterações máximo atingido.\n')
  
  return np.array(x_pontos), np.array(y_pontos), np.array(grad_f)

Testes¶

Nesta parte algumas funções exemplo serão testadas utilizando o Método de Descida do Gradiente com busca linear de Armijo.

Exemplo 1 - Função Esfera¶

O primeiro exemplo mais simples é com a Função Esfera, definida como:

$$ $$\begin{equation} f(\mathbf{x}) = \sum_{i=1}^{n} x^{2}_{i} \end{equation}$$ $$

Com $\mathbf{x} \in \mathbb{R}^{n}$, $n \in \mathbb{N}$. Foi considerado para os exemplos abaixo a função com 2 variáveis, ou seja, $n=2$.

In [7]:
f_esfera = lambda xk: xk[0]**2 + xk[1]**2
f_esfera_grad = nd.Gradient(f_esfera)
f_esfera_hessian = nd.Hessian(f_esfera)
plot_function(f_esfera, title = 'Função Esfera', ndim = 3)

Cujos resultados são os seguintes:

In [8]:
chute = [4, -3]
xi, yi, grad_norm = GradientDescent(f_esfera, f_esfera_grad, chute = chute, tol=1e-10)
df1_1 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||']) 
df1_1
Chute inicial: y = 25, x = [4, -3] 

Out[8]:
x1 x2 f(x1,x2) ||∇f||
0 4.00000000 -3.00000000 25.00000000 10.00000000
1 0.00000000 0.00000000 0.00000000 0.00000000
In [9]:
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados ***************

Número de iterações:  2

***** Primeira iteração *****

f([ 4 -3]) = 25.0

||∇f([ 4 -3])|| =  10.0


***** Última iteração: 2° *****

f([0. 0.]) = 0.0

||∇f([0. 0.])|| =  0.0

O algoritmo convergiu de acordo com critério de parada $\lVert\nabla f(\mathbf{x}) \rVert \leq 10^{-10}$, estabelecido anteriormente com $2$ iterações.

Abaixo temos 2 $\textit{plots}$, o primeiro com as curvas de nível e o caminho tomado pelo algoritmo durante a otimização, e o segundo com o valor da função custo a cada iteração.

In [10]:
plot_results(f_esfera, xi, yi, method = '1')

Exemplo 2 - Função Booth¶

A Função Booth é uma função de $2$ dimensões definida como:

$$ $$\begin{equation} f(\mathbf{x}) = (x_{1} + 2x_{2} - 7)^{2} + (2x_{1} + x_{2} - 5)^{2} \end{equation}$$ $$

Primeiro, temos seu gráfico:

In [11]:
f_booth = lambda xk: (xk[0] + 2*xk[1] - 7)**2 + (2*xk[0] + xk[1] - 5)**2
f_booth_grad = nd.Gradient(f_booth)
f_booth_hessian = nd.Hessian(f_booth)
plot_function(f_booth, title = 'Função Booth', ndim = 3)
In [12]:
chute = [5, 5]
xi, yi, grad_norm = GradientDescent(f_booth, f_booth_grad, chute = chute, tol=1e-10)
df2_1 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||']) 
Chute inicial: y = 164, x = [5, 5] 

In [13]:
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados ***************

Número de iterações:  61

***** Primeira iteração *****

f([5 5]) = 164.0

||∇f([5 5])|| =  76.4198927


***** Última iteração: 61° *****

f([1. 3.]) = 0.0

||∇f([1. 3.])|| =  0.0

Que atingiu o critério de parada após $61$ iterações do Algoritmo. Abaixo temos os resultados em gráficos:

In [14]:
plot_results(f_booth, xi, yi, dim = [np.linspace(-10, 15, 500), np.linspace(-10, 15, 500)], method = '1')
In [ ]:
 

Exemplo 3 - Função Styblinski-Tang com n = 2¶

Neste outro caso temos a Função Styblinski-Tang, definida como:

$$ $$\begin{equation} f(\mathbf{x}) = \dfrac{1}{2} \sum_{i = 1}^{n}(x_{i}^{4}-16x_{i}^{2}+5x_{i}) \end{equation}$$ $$

Que foi considerada com $n = 2$, ou seja, duas dimensões e avaliada em $x_{0} = [-3, -4]$

In [15]:
f_tang = lambda xk: (xk[0]**4 - 16*xk[0]**2 + 5*xk[0] + xk[1]**4 - 16*xk[1]**2 + 5*xk[1])/2
f_tang_grad = nd.Gradient(f_tang)
f_tang_hessian = nd.Hessian(f_tang)
plot_function(f_tang, title = 'Função Styblinski-Tang com n = 2', ndim = 3, angle = (40, 50))
In [16]:
chute = [-3, -4]
xi, yi, grad_norm = GradientDescent(f_tang, f_tang_grad, chute = chute, tol=1e-7)
df3_1 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||']) 
df3_1
Chute inicial: y = -49.0, x = [-3, -4] 

Out[16]:
x1 x2 f(x1,x2) ||∇f||
0 -3.00000000 -4.00000000 -49.00000000 61.59951299
1 -2.89062500 -2.07812500 -69.58226913 17.80632896
2 -2.91834593 -3.19067526 -76.76195048 11.42561239
3 -2.90221876 -2.83398927 -78.25061301 2.32193091
4 -2.90363926 -2.90653571 -78.33217526 0.10402795
5 -2.90352553 -2.90328682 -78.33233035 0.00855307
6 -2.90353471 -2.90355395 -78.33233140 0.00068935
7 -2.90353397 -2.90353242 -78.33233141 0.00005565
8 -2.90353403 -2.90353416 -78.33233141 0.00000449
9 -2.90353403 -2.90353402 -78.33233141 0.00000036
10 -2.90353403 -2.90353404 -78.33233141 0.00000042
11 -2.90353403 -2.90353399 -78.33233141 0.00000140
12 -2.90353403 -2.90353403 -78.33233141 0.00000011
13 -2.90353403 -2.90353402 -78.33233141 0.00000013
14 -2.90353403 -2.90353403 -78.33233141 0.00000015
15 -2.90353403 -2.90353402 -78.33233141 0.00000018
16 -2.90353403 -2.90353403 -78.33233141 0.00000008

Que atingiu o critério de tolerância $\leq 10^{-7}$ após $11$ iterações do Algoritmo.

In [17]:
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados ***************

Número de iterações:  17

***** Primeira iteração *****

f([-3 -4]) = -49.0

||∇f([-3 -4])|| =  61.59951299


***** Última iteração: 17° *****

f([-2.90353403 -2.90353403]) = -78.33233141

||∇f([-2.90353403 -2.90353403])|| =  8e-08
In [18]:
plot_results(f_tang, xi, yi, dim = [np.linspace(-5, 5, 500), np.linspace(-5, 5, 500)], method = '1')
In [ ]:
 

Exemplo 4 - Função Rosenbrock com n = 2¶

Outro exemplo é a Função Rosenbrock, que é uma função $n-dimensional$ unimodal (com apenas 1 mínimo) sendo que esse mínimo se encontra em um vale parabólico estreito $[2]$.

$$ $$

A função é definida como:

$$ $$\begin{equation} f(\mathbf{x}) = \sum_{i=1}^{n-1} \big(100 (x_{i+1} - x_ {i}^{2})^{2} + (x_{i}-1)^{2}\big) \end{equation}$$ $$

Foi considerada a função com $n = 2$ e a constante (100) que multiplica o primeiro termo da Função Rosenbrock foi removido para facilitar a visualização com valores menores.

$$ $$

Os resultados abaixo referem-se ao Método do Gradiente calculado com chute inicial $x_{0} = [-5, -25]$.

In [19]:
f_rosenbrock = lambda xk: (xk[1] - xk[0]**2)**2 + (1 - xk[0])**2
f_rosenbrock_grad = nd.Gradient(f_rosenbrock)
plot_function(f_rosenbrock, title = 'Função Rosenbrock $n=2$', ndim = 3, dom = np.linspace(-2, 2, 500), angle = (30, 40))
In [20]:
chute = [-5, -25]
xi, yi, grad_norm = GradientDescent(f_rosenbrock, f_rosenbrock_grad, chute = chute, tol=1e-10)
df4_1 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Chute inicial: y = 2536, x = [-5, -25] 


*************** Resultados ***************

Número de iterações:  335

***** Primeira iteração *****

f([ -5 -25]) = 2536.0

||∇f([ -5 -25])|| =  1016.9287094


***** Última iteração: 335° *****

f([1. 1.]) = 0.0

||∇f([1. 1.])|| =  0.0

Atingindo o critério desejado após 335 iterações.

In [21]:
plot_results(f_rosenbrock, xi, yi, dim = [np.linspace(-13, 4, 1000), np.linspace(-200, 250, 1000)], method = '1')

Podemos perceber que o método oscilou muito até que fosse atingido o critério de parada pré-estabelecido. Variando o critério de tolerância para uma precisão muito maior, pode-se perceber que o método não converge, ou seja, essas oscilações permanecem e não atingem o desejado. Esse resultado será importante para se comparar com os outros Algoritmos como o Método de Newton e o BFGS, que serão apresentados abaixo e que são teoricamente melhores que o Método de Descida do Gradiente.

In [ ]:
 

Exemplo 5 - Função McCormick¶

A Função McCormick, por sua vez, é uma função de 2 variáveis reais e também unimodal, definida como:

$$ $$\begin{equation} f(\mathbf{x}) = \sin{(x_{1}+x_{2})} + (x_{1}-x_{2})^{2} - 1.5x_{1} + 2.5x_{2} + 1 \end{equation}$$ $$

Essa função foi avaliada com $x_{0} = [-7, 6]$.

In [22]:
f_cormick = lambda xk: np.sin(xk[0] + xk[1]) + (xk[0] - xk[1])**2 - 1.5*xk[0] + 2.5*xk[1] + 1
f_cormick_grad = nd.Gradient(f_cormick)
plot_function(f_cormick, title = 'Função McCormick', ndim = 3, dom = np.linspace(-4, 4, 500), angle = (30, 50))
In [23]:
chute = [-7, 6]
xi, yi, grad_norm = GradientDescent(f_cormick, f_cormick_grad, chute = [-7, 6], tol=1e-10)  
df5_1 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Chute inicial: y = 194.6585290151921, x = [-7, 6] 


*************** Resultados ***************

Número de iterações:  45

***** Primeira iteração *****

f([-7  6]) = 194.65852902

||∇f([-7  6])|| =  39.62530073


***** Última iteração: 45° *****

f([-0.54719755 -1.54719755]) = -1.91322295

||∇f([-0.54719755 -1.54719755])|| =  0.0

Que satisfez o critério de parada com $\lVert \nabla f(x_{k}) \rVert \leq 10^{-10}$ após $46$ iterações.

In [24]:
plot_results(f_cormick, xi, yi, dim = [np.linspace(-45, 8, 1000), np.linspace(-45, 8, 1000)], method = '1')
In [ ]:
 

Método de Newton com Busca de Armijo¶

O Método de Newton, diferentemente do método apresentado acima, consiste em se obter uma direção $d_k$ tal que:

\begin{equation} d_k = - (\nabla^{2}f(\vec{x}))^{-1} \cdot \nabla f(\vec{x}) \end{equation}$$ $$

Em que $\nabla^{2}f(\vec{x}) = \mathbb{H}(x)_{(ij)} \equiv \frac{\partial^{2} f}{\partial x_{i} \partial x_{j} }$ é a matriz Hessiana.

Um fato a se destacar é que longe do ponto, $\nabla^{2}f(x_{k})$ pode ser não definida positiva. Caso isso ocorra, a direção do Método do Gradiente deve ser empregado como direção de descida.

A partir de então, calcular o novo ponto $x_{k+1}$ em que a norma do gradiente $\lVert \nabla f(x_{k+1}) \rVert$ seja menor e o custo da função $f(x_{k+1})$ diminua. O algoritmo esquematizado está abaixo:



Algoritmo 3 - Método de Newton com condição de Armijo




Dado um ponto $x_{0} \in \Re^{n}$ ($\textit{initial guess}$), para $\nabla f(x_{k}) \neq 0$ e $\nabla^{2}f(x_{k})$ matriz não singular, realizar os seguintes passos:

1. Calcular a direção de descida $d_{k}$ resolvendo o sistema $\nabla^{2}f(x_{k}) \, d_{k} = - \nabla f(x_{k})$;

2. Verificar se $d_{k}^{T}\nabla f(x_{k}) \leq -0.001 \lVert d_{k} \rVert \lVert \nabla f(x_{k}) \rVert$. Se não for satisfeita, usar $d_{k} = - \nabla f(x_{k})$ (Método do Gradiente);

3. Realizar Busca de Armijo e encontrar um $\lambda_k$ que satisfaça a Desigualdade 1;

4. Calcular o novo ponto $x_{k+1} = x_{k} + \lambda_{k} d_{k}$;




Em que tais etapas são realizadas com os mesmos critérios de parada dados no Método de Descida do Gradiente.

Ademais, a convergência do Método de Newton é quadrática, ou seja, a sequência {$x_k$} converge para $x^{*}$ quando:

\begin{equation} \lVert x_{k+1} - x^{*} \rVert \leq \epsilon \, \lVert x_{k} - x^{*} \rVert^{2} \end{equation}

$\forall k \in \mathbb{N}$ e com $\epsilon > 0$.

Sendo assim teoricamente melhor que o método acima e o método apresentado a seguir.

A chamada para a Função que implementa o Método de Newton com Condição de Armijo é:

Newton(f, f_grad, f_hessian, chute, sigma, tol)

Onde:

$$ $$
  • f é a função desejada (tipo: callable)

  • f_grad é o gradiente da função f (tipo: callable)

  • f_hessian é a hessiana da função f (tipo: callable)

  • chute é o $x_{0}$, ponto inicial (tipo: np.array ou list)

  • sigma é a constante de decréscimo de Armijo (tipo: float $\in (0, 1)$. Padrão $= 0.02$)

  • tol é a tolerância desejada ($\lVert \nabla f(x_{k+1}) \rVert \leq 10^{-8}$) por exemplo (tipo: float $> 0$. Padrão = 1e-8 = $10^{-8}$)

$$ $$

E pode retornar até 3 argumentos que são tipo np.array utilizados para fazer gráficos e outras análises visuais.

$$ $$
  • x são pontos $x_k$ calculados pelo Método

  • y são valores de $f(x_{k})$ calculados pelo Método

  • grad_f são valores de $\lVert \nabla f(x_k) \rVert$ calculados pelo Método

$$ $$

Para suprimir algum retorno, por exemplo, o grad_f, basta fazer:

$$ $$
x, y, _ = Newton(f, f_grad, f_hessian, chute, sigma, tol)
$$ $$

Exemplos de como utilizar tais funções serão fornecidos mais adiante com funções pré-estabelecidas. Abaixo temos a implementação do Algoritmo 3:

Algoritmo¶

In [25]:
def Newton(f, f_grad, f_hessian, chute, sigma = 0.02, tol=1e-8):

  """Algoritmo de Descida do Gradiente com busca de Armijo. Esquema:

        I - Direção de descida: d_k := −∇f(x).
        II - Determinação do passo com busca de Armijo.
        III - Obtém o próximo candidato.

        Parâmetros
        --------------------
        f : callable
            Função custo
        f_grad : callable
            Gradiente da função f
        f_hessian : callable
            Hessiana da função f 
        chute : array
            Valor inicial de x ("chute")
        sigma : float, opcional
            Constante de Decréscimo de Armijo. Padrão = 0.02
        tol : float, opcional
            Tolerância desejada. Padrão = 1e-8
        
        Saída
        --------------------
        xk : array
            valores de xk
        yk : array
            valores da função custo e cada xk
        grad_f : array
            valores da norma do gradiente da f em xk"""

  # Valores iniciais de xk, fk, grad_fk e max_iter
  xk = chute    
  fk = f(xk)
  grad_fk = f_grad(xk)
  grad_fk_norm = np.linalg.norm(grad_fk)
  max_iter = 500

  # Inicializa o número de iterações e a lista para fazer os plots dos valores de x e y
  num_iter = 0
  x_pontos = [xk]
  y_pontos = [fk]
  grad_f = [grad_fk_norm]

  f_hessian_xk = f_hessian(xk)  
  
  # Calcula nova iteração com busca de Armijo
  while (grad_fk_norm > tol and num_iter < max_iter):
        
    f_hessian_xk = f_hessian(xk)

    if (np.linalg.norm(f(xk)) > 1e16) or (grad_fk_norm > 1e16)  or  (np.linalg.norm(f_hessian_xk) > 1e16):
      print("\nErro: Overflow\n")
      break 
    
    # Não invertível
    if abs(np.linalg.det(f_hessian_xk)) <= 1e-3:
      break 

    # Determina a direção resolvendo o sistema:  \nabla^{2} f(x_{k}) . p_{k} = - \nabla f(x_{k})
    pk = -np.linalg.solve(f_hessian_xk, grad_fk)
    lambda_ = 1 

    # Se estiver longe do ponto, a matriz Hessiana pode ser não definida positiva.
    # Se isso acontecer, seguir pelo Método do Gradiente:
    if not (np.dot(pk, grad_fk) < (-0.001 * np.linalg.norm(f_grad(xk)) * np.linalg.norm(pk))):
      
      # Descida do gradiente
      pk = -f_grad(xk)
      
      # Faz a busca de Armijo, obtendo o passo lambda e a função custo naquele passo
      lambda_, fk = Armijo_Search(f, xk, pk, sigma=sigma)

    else:
      lambda_, fk = Armijo_Search(f, xk, pk, sigma=sigma)
    
    # calcula x_{k+1}
    xk = xk + lambda_ * pk
    grad_fk = f_grad(xk)
    grad_fk_norm = np.linalg.norm(grad_fk)

    # Itera mais uma vez 
    num_iter += 1
    x_pontos.append(xk)
    y_pontos.append(f(xk))
    grad_f.append(grad_fk_norm)

    # número máximo de iterações
    if num_iter == max_iter:
      print('\nNúmero máximo de iterações atingido.\n')
      break

  if abs(np.linalg.det(f_hessian_xk)) <= 1e-3:
    print("\nMatriz Hessiana não invertível. Mudar o ponto inicial.\n") 
    return [], [], []
  
  elif np.min(np.linalg.eigvals(f_hessian_xk)) <= 1e-3:
    print("\nMatriz Hessiana não PD. Mudar o ponto inicial. Ponto de sela encontrado.\n")
  
  
  return np.array(x_pontos), np.array(y_pontos), np.array(grad_f)

Testes¶

Nesta parte algumas funções exemplo serão testadas utilizando o Método de Newton com busca linear de Armijo, como os exemplos 3, 4 e 5 do Método de Descida do Gradiente.

Exemplo 1 - Função Esfera¶

Testando o Método de Newton na Função Esfera definida anteriormente, temos:

In [26]:
f_esfera = lambda xk: xk[0]**2 + xk[1]**2
f_esfera_grad = nd.Gradient(f_esfera)
f_esfera_hessian = nd.Hessian(f_esfera)
plot_function(f_esfera, title = 'Função Esfera', ndim = 3)
In [27]:
chute = [4, -3]
xi, yi, grad_norm =  Newton(f_esfera, f_esfera_grad, f_esfera_hessian, chute = chute, tol=1e-10)
df1_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||']) 
df1_2
Out[27]:
x1 x2 f(x1,x2) ||∇f||
0 4.00000000 -3.00000000 25.00000000 10.00000000
1 -0.00000000 0.00000000 0.00000000 0.00000000
In [28]:
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados ***************

Número de iterações:  2

***** Primeira iteração *****

f([ 4 -3]) = 25.0

||∇f([ 4 -3])|| =  10.0


***** Última iteração: 2° *****

f([-0.  0.]) = 0.0

||∇f([-0.  0.])|| =  0.0

O Método de Newton convergiu após 2 iterações. Abaixo temos os $\textit{plots}$ com mais alguns resultados:

In [29]:
plot_results(f_esfera, xi, yi, method = '2')
In [ ]:
 

Exemplo 2 - Função Booth¶

Para a Função Booth temos os seguintes resultados:

In [30]:
f_booth = lambda xk: (xk[0] + 2*xk[1] - 7)**2 + (2*xk[0] + xk[1] - 5)**2
f_booth_grad = nd.Gradient(f_booth)
f_booth_hessian = nd.Hessian(f_booth)
plot_function(f_booth, title = 'Função Booth', ndim = 3)
In [31]:
chute = [5, 5]
xi, yi, grad_norm = Newton(f_booth, f_booth_grad, f_booth_hessian, chute = chute, tol=1e-10)
df2_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||']) 
df2_2
Out[31]:
x1 x2 f(x1,x2) ||∇f||
0 5.00000000 5.00000000 164.00000000 76.41989270
1 1.00000000 3.00000000 0.00000000 0.00000000
In [32]:
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados ***************

Número de iterações:  2

***** Primeira iteração *****

f([5 5]) = 164.0

||∇f([5 5])|| =  76.4198927


***** Última iteração: 2° *****

f([1. 3.]) = 0.0

||∇f([1. 3.])|| =  0.0

Que atingiu o critério de parada $\lVert\nabla f(\mathbf{x}) \rVert \leq 10^{-10}$ após 2 iterações. Aqui já vemos uma diferença considerável se comparada com o Método de Descida do Gradiente, uma vez que a convergência só foi atingida após $61$ iterações, enquanto que no Método de Newton foi praticamente imediata sob as mesmas condições de partida e parada.

In [33]:
plot_results(f_booth, xi, yi, dim = [np.linspace(-10, 15, 500), np.linspace(-10, 15, 500)], method = '2')

Exemplo 3 - Função Styblinski-Tang com n = 2¶

Neste caso temos novamente a Função Styblinski-Tang com $n = 2$ com um caso bastante interessante do Método de Newton. Primeiro, o algoritmo será executado com um passo inicial (chute inicial) de $x_{0} = [0, 2]$.

In [34]:
f_tang = lambda xk: (xk[0]**4 - 16*xk[0]**2 + 5*xk[0] + xk[1]**4 - 16*xk[1]**2 + 5*xk[1])/2
f_tang_grad = nd.Gradient(f_tang)
f_tang_hessian = nd.Hessian(f_tang)
plot_function(f_tang, title = 'Função Styblinski-Tang com n = 2', ndim = 3, angle = (40, 50))
In [35]:
chute = [0, 2]
xi, yi, grad_norm = Newton(f_tang, f_tang_grad, f_tang_hessian, chute = chute, tol=1e-7)
df3_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Matriz Hessiana não PD. Mudar o ponto inicial. Ponto de sela encontrado.


*************** Resultados ***************

Número de iterações:  5

***** Primeira iteração *****

f([0 2]) = -19.0

||∇f([0 2])|| =  13.72953022


***** Última iteração: 5° *****

f([0.15673126 2.74680277]) = -24.8338343

||∇f([0.15673126 2.74680277])|| =  0.0

O Algoritmo 3 identificou que a condição sobre a Matriz Hessiana ser Positiva Definida foi violada. Trata-se de um ponto de sela.

In [36]:
plot_results(f_tang, xi, yi, dim = [np.linspace(-5, 5, 500), np.linspace(-5, 5, 500)], method = '2')

Esse é um clássico problema do Método de Newton, tanto com $1$ variável real como multivariado, a escolha incial é importante para a convergência do método.

$$ $$

Mudando o ponto de começo $x_{0}$ para $x_{0} = [-3, -4]$, temos os seguintes resultados:

In [37]:
chute = [-3, -4]
xi, yi, grad_norm = Newton(f_tang, f_tang_grad, f_tang_hessian, chute = [-3, -4], tol=1e-7)
df3_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||']) 
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados ***************

Número de iterações:  6

***** Primeira iteração *****

f([-3 -4]) = -49.0

||∇f([-3 -4])|| =  61.59951299


***** Última iteração: 6° *****

f([-2.90353403 -2.90353403]) = -78.33233141

||∇f([-2.90353403 -2.90353403])|| =  0.0

Chegando no mesmo valor que o Método do Gradiente com apenas 6 passos.

In [38]:
plot_results(f_tang, xi, yi, dim = [np.linspace(-5, 5, 500), np.linspace(-5, 5, 500)], method = '2')

Que convergiu para o mínimo da função.

In [ ]:
 

Exemplo 4 - Função Rosenbrock com n = 2¶

A função testada é novamente a Função Rosenbrock com n = 2 utilizada no Exemplo do Método anterior.

Os resultados abaixo referem-se ao Método de Newton calculado com chute inicial $x_{0} = [-5, -25]$.

In [39]:
f_rosenbrock = lambda xk: (xk[1] - xk[0]**2)**2 + (1 - xk[0])**2
f_rosenbrock_grad = nd.Gradient(f_rosenbrock)
f_rosenbrock_hessian = nd.Hessian(f_rosenbrock)
plot_function(f_rosenbrock, title = 'Função Rosenbrock $n=2$', ndim = 3, dom = np.linspace(-2, 2, 500), angle = (30, 40))
In [40]:
chute = [-5, -25]
xi, yi, grad_norm = Newton(f_rosenbrock, f_rosenbrock_grad, f_rosenbrock_hessian, chute = chute, tol=1e-10) 
df4_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados ***************

Número de iterações:  13

***** Primeira iteração *****

f([ -5 -25]) = 2536.0

||∇f([ -5 -25])|| =  1016.9287094


***** Última iteração: 13° *****

f([1. 1.]) = 0.0

||∇f([1. 1.])|| =  0.0

No Método do Gradiente, para esse problema, o critério de parada $\lVert\nabla f(\mathbf{x}) \rVert \leq 10^{-10}$ foi satisfeito somente após $335$ iterações. O Método de Newton, por sua vez, sob as mesmas condições, atingiu o critério após apenas $13$ iterações, sendo outra diferença notável da convergência dos dois métodos.

In [41]:
plot_results(f_rosenbrock, xi, yi, dim = [np.linspace(-13, 4, 1000), np.linspace(-200, 250, 1000)], method = '2')

Além disso, o Método de Newton não oscilou como o Método do Gradiente para encontrar a solução em que o critério fosse satisfeito, que acaba sendo a solução da Função Rosenbrock bi-dimensional.

In [ ]:
 

Exemplo 5 - Função McCormick¶

A Função McCormick foi avaliada com $x_{0} = [-7, 6]$ e com critério $\lVert\nabla f(\mathbf{x}) \rVert \leq 10^{-10}$, como anteriormente.

In [42]:
f_cormick = lambda xk: np.sin(xk[0] + xk[1]) + (xk[0] - xk[1])**2 - 1.5*xk[0] + 2.5*xk[1] + 1
f_cormick_grad = nd.Gradient(f_cormick)
f_cormick_hessian = nd.Hessian(f_cormick)
plot_function(f_cormick, title = 'Função McCormick', ndim = 3, dom = np.linspace(-4, 4, 500), angle = (30, 50))
In [43]:
chute = [-7, 6]
xi, yi, grad_norm = Newton(f_cormick, f_cormick_grad, f_cormick_hessian, chute = chute, tol=1e-10)  
df5_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados ***************

Número de iterações:  5

***** Primeira iteração *****

f([-7  6]) = 194.65852902

||∇f([-7  6])|| =  39.62530073


***** Última iteração: 5° *****

f([-0.54719755 -1.54719755]) = -1.91322295

||∇f([-0.54719755 -1.54719755])|| =  0.0

Novamente encontrando os mesmos valores do Método de Newton em muito menos iterações do Algoritmo.

In [44]:
plot_results(f_cormick, xi, yi, dim = [np.linspace(-45, 8, 1000), np.linspace(-45, 8, 1000)], method = '2')

Exemplo 6 - Função Colville¶

A Função Colville é uma função de $4$ dimensões, definida como:

$$ $$\begin{equation} f(\mathbf x) = 100(x_{1}^2-x_{2})^2 + (x_{1}-1)^2 + (x_{3} -1)^2 + 90(x_{3}^2-x_{4})^2\\+10.1((x_{2} - 1)^2+ (x_{4}-1)^2) +19.8(x_{2} -1)(x_{4} -1) \end{equation}$$ $$

Que será avaliada no ponto $x_{0} = [3, 1, 5, 4]$:

In [45]:
def f_colville(xk):

  """
    Função Colville - 4 dimensões: https://www.sfu.ca/~ssurjano/colville.html

     Parâmetros
    --------------------
    xk : np.array or list
        array de pontos xk
    
    Saída
    --------------------
    y : array
        valores da função custo
  """

  x1 = xk[0]
  x2 = xk[1]
  x3 = xk[2]
  x4 = xk[3]

  term1 = 100 * (x1**2-x2)**2;
  term2 = (x1-1)**2;
  term3 = (x3-1)**2;
  term4 = 90 * (x3**2-x4)**2;
  term5 = 10.1 * ((x2-1)**2 + (x4-1)**2);
  term6 = 19.8*(x2-1)*(x4-1);

  y = term1 + term2 + term3 + term4 + term5 + term6

  return y
In [46]:
f_colville_grad = nd.Gradient(f_colville)
f_colville_hessian = nd.Hessian(f_colville)
In [47]:
chute = [3, 1, 5, 4]
xi, yi, grad_norm = Newton(f_colville, f_colville_grad, f_colville_hessian, chute = chute, tol=1e-8)  
df6_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'x3':xi[:,2], 'x4':xi[:,3], 'f(x1,x2,x3,x4)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'x3', 'x4', 'f(x1,x2,x3,x4)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados ***************

Número de iterações:  15

***** Primeira iteração *****

f([3 1 5 4]) = 46200.9

||∇f([3 1 5 4])|| =  39215.92871169


***** Última iteração: 15° *****

f([1. 1. 1. 1.]) = 0.0

||∇f([1. 1. 1. 1.])|| =  0.0
In [48]:
df6_2
Out[48]:
x1 x2 x3 x4 f(x1,x2,x3,x4) ||∇f||
0 3.00000000 1.00000000 5.00000000 4.00000000 46200.90000000 39215.92871169
1 2.29966013 2.93255387 4.49767692 17.87089851 4627.45416071 4399.11465877
2 1.66303161 1.70962687 3.09191394 6.84928366 1210.45916989 3123.64881103
3 1.02829368 0.25151366 2.56382101 5.83413244 286.58824030 767.13016041
4 0.82740311 0.48596161 1.92643823 3.12318333 62.53754579 421.28883716
5 0.54946097 0.15520747 1.64994513 2.56335455 10.79500134 102.82494463
6 0.50791572 0.24068237 1.41680691 1.93193006 1.54159644 40.57739563
7 0.58543341 0.33528193 1.30739964 1.69270586 0.48876832 8.82368418
8 0.79961196 0.59452106 1.19068036 1.40342568 0.36160394 18.08480041
9 0.83862492 0.70294929 1.13752332 1.29113716 0.08062932 1.56255741
10 0.99250921 0.96141233 1.02560337 1.03918358 0.07177742 11.74202871
11 0.99147987 0.98296381 1.00919960 1.01806179 0.00030743 0.17871781
12 0.99982048 0.99957110 1.00026150 1.00044089 0.00000127 0.04548513
13 0.99999827 0.99999652 1.00000180 1.00000352 0.00000000 0.00003561
14 1.00000000 1.00000000 1.00000000 1.00000000 0.00000000 0.00000000

Convergindo para o mínimo global $f(1,1,1,1) = 0$ após $15$ passos.

Método Broyden–Fletcher–Goldfarb–Shanno (BFGS) com condição de Armijo¶

O método $\textbf{B}$royden–$\textbf{F}$letcher–$\textbf{G}$oldfarb–$\textbf{S}$hanno ($\textbf{BFGS}$) é um método Quasi-Newton, sendo uma alternativa apresentada para o problema encontrado no Método de Newton apresentado acima, em que é necessário calcular a matriz Hessiana $\nabla^{2}{f(\vec{x})}$, ($x \in \Re^{n}$).

Tal problema é contornado pelo cálculo de uma matriz que serve como uma aproximação para a matriz Hessiana, não sendo necessários, portanto, o seu cálculo diretamente e o cálculo de sua inversa, como foi realizado no método acima.

O algoritmo e sua esquematização foram inspiradas em $[4]$.



Algoritmo 4 - Método BFGS com condição de Armijo




A partir de um ponto $x_{0} \in \Re^{n}$ ($\textit{initial guess}$) e uma matriz inicial $B_{0} \equiv \mathbb{I}_{n}$, com $n \in \mathbb{N}$, realizar os seguintes passos:

$$ $$

1. Obter uma direção de descida $d_{k} = - B_{k} \cdot \nabla f(x_{k})$;

$$ $$

2. Realizar busca linear e determinar $\lambda_{k}$ que satisfaz a condição de Armijo (Desigualdade 1);

$$ $$

3. Fazer $s_{k} = \lambda_{k} \, d_{k}$ e atualizar o ponto $x_{k+1} = x_{k} + \lambda_{k} \, d_{k} = x_{k} + s_{k}$;

$$ $$

4. Calcular a diferença de gradientes $y_{k} = \nabla f(x_{k+1}) - \nabla f(x_{k})$;

$$ $$

5. Atualizar a matriz $B_{k+1}$ que aproxima a Hessiana:

\begin{equation}\tag{2} B_{k+1} = B_{k} + \dfrac{(s^{T}_{k} y_{k}+ y^{T}_{k}B_{k}y_{k})(s_{k}s^{T}_{k})}{(s^{T}_{k}y_{k})^{2}} - \dfrac{B_{k}y_{k}s^{T}_{k} + s_{k}y^{T}_{k}B_{k}}{s^{T}_{k}y_{k}} \end{equation}$$ $$



$$ $$

Em que tais etapas são realizadas até os critérios apresentados no começo deste trabalho forem satisfeitos.

$$ $$

Vale ressaltar ainda que a convergência do Método BFGS é superlinear, ou seja, a sequência {$x_n$} converge para $x^{*}$ quando:

\begin{equation} \lim_{n→\infty} \dfrac{\lVert x_{n+1} - x^{*} \rVert}{\lVert x_{n} - x^{*} \rVert} = 0 \end{equation}$$ $$

Sendo assim, teoricamente melhor que o Método de Descida do Gradiente, porém pior que o Método de Newton.

A chamada para a Função que implementa o Método BFGS com Condição de Armijo é:

BFGS(f, f_grad, chute, sigma, tol)

Onde:

$$ $$
  • f é a função desejada (tipo: callable)

  • f_grad é o gradiente da função f (tipo: callable)

  • chute é o $x_{0}$, ponto inicial (tipo: np.array ou list)

  • sigma é a constante de decréscimo de Armijo (tipo: float $\in (0, 1)$. Padrão $= 0.02$)

  • tol é a tolerância desejada ($\lVert \nabla f(x_{k+1}) \rVert \leq 10^{-8}$) por exemplo (tipo: float $> 0$. Padrão = $10^{-8}$)

$$ $$

E pode retornar até 3 argumentos que são tipo np.array utilizados para fazer gráficos e outras análises visuais.

$$ $$
  • x são pontos $x_k$ calculados pelo Método

  • y são valores de $f(x_{k})$ calculados pelo Método

  • grad_f são valores de $\lVert \nabla f(x_k) \rVert$ calculados pelo Método

$$ $$

Para suprimir algum retorno, por exemplo, o grad_f, basta fazer:

$$ $$
x, y, _ = BFGS(f, f_grad, chute, sigma, tol)
$$ $$

Exemplos de como utilizar tais funções serão fornecidos mais adiante com funções pré-estabelecidas. Abaixo temos a implementação do Algoritmo 4:

Algoritmo¶

In [49]:
def BFGS(f, f_grad, chute, sigma = 0.02, tol=1e-8):

  """ Algoritmo Broyden–Fletcher–Goldfarb–Shanno (BFGS) com busca de Armijo.

    Parâmetros
    --------------------
    f : callable
        Função custo
    f_grad : callable
        Gradiente da função f
    chute : array
        Valor inicial de x ("chute")
    sigma : float, opcional 
        Constante de decréscimo de Armijo. Padrão = 0.02
    tol : float, opcional
        Tolerância desejada. Padrão = 1e-8
        
    Saída
    --------------------
    xk : array
        valores de xk
    yk : array
        valores da função custo em cada xk  
    grad_f : array
      valores da norma do gradiente da f em xk"""


  # dimensão do problema e matriz B_{k = 0} = Id(n) (identidade de ordem n)
  n = len(chute)
  Bk =  np.eye(len(chute))
  max_iter = 500

  # Valores iniciais de x_{k+1}, x_{k}, f_{x_{k}}, grad_{f(x_{k})} e norma do gradiente
  xk = chute    
  xprev = chute
  fk = f(xk)
  grad_fk = f_grad(xk)
  grad_fk_norm = np.linalg.norm(grad_fk)

  # Inicializa o número de iterações e a lista para fazer os plots dos valores de x e y
  num_iter = 0
  x_pontos = [xk]
  y_pontos = [fk]
  grad_f = [grad_fk_norm]

  print(f'Chute inicial: y = {fk}, x = {xk} \n')

  while (grad_fk_norm > tol and num_iter < max_iter):
    
    if (np.linalg.norm(f(xk)) > 1e16) or (grad_fk_norm > 1e16):
      print("\nErro: Overflow\n")
      break 

    # direção (método de Newton)
    dk = - Bk @ f_grad(xprev) 

    # performa busca de Armijo com x_{k} (x prévio)
    lambda_, _ = Armijo_Search(f, xprev, dk, sigma)
    
    # direção e valor de x_{k+1} = x_{k} + \lambda * d_{k} = x_{k} + s_{k} 
    sk = lambda_ * dk
    xk = xk + lambda_ * dk

    # y_{k} é a diferença do gradiente da f em x_{k+1} com o gradiente da f em x_{k} 
    yk = np.subtract(f_grad(xk), f_grad(xprev)) 
    
    # arruma os tamanhos 
    yk = np.array([yk]).reshape(-1, 1)
    sk = np.array([sk]).reshape(-1, 1)

    # denominador do BFGS update
    denominador = np.dot(sk.T, yk)[0]  # escalar
    
    # numerador da segunda parcela do BFGS update
    num_1 = ((denominador + (yk.T @ (Bk @ yk))) * (sk @ sk.T))/(denominador ** 2) 

    # numerador da terceira parcela do BFGS update
    num_2 = ((Bk @ (yk @ sk.T)) + (sk @ (yk.T @ Bk)))/denominador

    # BFGS update
    Bk = Bk + np.subtract(num_1, num_2)

    grad_fk = f_grad(xk)
    grad_fk_norm = np.linalg.norm(grad_fk)

    x_pontos.append(xk)
    y_pontos.append(np.round(f(xk), 10))
    grad_f.append(grad_fk_norm)

    xprev = xk
    num_iter += 1

  # número máximo de iterações
  if num_iter == max_iter:
    print('\nNúmero máximo de iterações atingido.\n')

  return np.array(x_pontos), np.array(y_pontos), np.array(grad_f)

Testes¶

Nesta parte algumas funções exemplo serão testadas utilizando o Método BFGS com busca linear de Armijo.

Exemplo 1 - Função Styblinski-Tang com n = 2¶

In [50]:
f_tang = lambda xk: (xk[0]**4 - 16*xk[0]**2 + 5*xk[0] + xk[1]**4 - 16*xk[1]**2 + 5*xk[1])/2
f_tang_grad = nd.Gradient(f_tang)
f_tang_hessian = nd.Hessian(f_tang)
plot_function(f_tang, title = 'Função Styblinski-Tang com n = 2', ndim = 3, angle = (40, 50))
In [51]:
chute = [-3, -2]
xi, yi, grad_norm = BFGS(f_tang, f_tang_grad, chute = chute, tol=1e-7)
df1_3 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Chute inicial: y = -68.0, x = [-3, -2] 


*************** Resultados ***************

Número de iterações:  11

***** Primeira iteração *****

f([-3 -2]) = -68.0

||∇f([-3 -2])|| =  18.82817038


***** Última iteração: 11° *****

f([-2.90353403 -2.90353403]) = -78.33233141

||∇f([-2.90353403 -2.90353403])|| =  6e-08
In [52]:
plot_results(f_tang, xi, yi, dim = [np.linspace(-5, 5, 500), np.linspace(-5, 5, 500)], method = '3')

Nessa primeira função é interessante perceber que o número de iterações para o critério de parada foi maior que no Método de Newton.

Exemplo 2 - Função Rosenbrock com n = 2¶

Testando a função Rosenbrock $2$-$d$ no método BFGS, temos:

In [53]:
f_rosenbrock = lambda xk: (xk[1] - xk[0]**2)**2 + (1 - xk[0])**2
f_rosenbrock_grad = nd.Gradient(f_rosenbrock)
f_rosenbrock_hessian = nd.Hessian(f_rosenbrock)
plot_function(f_rosenbrock, title = 'Função Rosenbrock $n=2$', ndim = 3, dom = np.linspace(-2, 2, 500), angle = (30, 40))
In [54]:
chute = [-5, -25]
xi, yi, grad_norm = BFGS(f_rosenbrock, f_rosenbrock_grad, chute = [-5, -25], tol=1e-10)
df2_3 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Chute inicial: y = 2536, x = [-5, -25] 


*************** Resultados ***************

Número de iterações:  17

***** Primeira iteração *****

f([ -5 -25]) = 2536.0

||∇f([ -5 -25])|| =  1016.9287094


***** Última iteração: 17° *****

f([1. 1.]) = 0.0

||∇f([1. 1.])|| =  0.0

Percebe-se que o Método BFGS teve mais iterações que o Método de Newton e menos iterações que o Método do Gradiente, até que o critério de parada fosse atingido.

In [55]:
plot_results(f_rosenbrock, xi, yi, dim = [np.linspace(-13, 15, 1000), np.linspace(-200, 250, 1000)], method = '3')

Exemplo 3 - Função McCormick¶

In [56]:
f_cormick = lambda xk: np.sin(xk[0] + xk[1]) + (xk[0] - xk[1])**2 - 1.5*xk[0] + 2.5*xk[1] + 1
f_cormick_grad = nd.Gradient(f_cormick)
f_cormick_hessian = nd.Hessian(f_cormick)
plot_function(f_cormick, title = 'Função McCormick', ndim = 3, dom = np.linspace(-4, 4, 500), angle = (30, 50))
In [57]:
chute = [-7, 6]
xi, yi, grad_norm = BFGS(f_cormick, f_cormick_grad, chute = chute, tol=1e-10)  
df3_3 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)  
Chute inicial: y = 194.6585290151921, x = [-7, 6] 


*************** Resultados ***************

Número de iterações:  9

***** Primeira iteração *****

f([-7  6]) = 194.65852902

||∇f([-7  6])|| =  39.62530073


***** Última iteração: 9° *****

f([-0.54719755 -1.54719755]) = -1.91322296

||∇f([-0.54719755 -1.54719755])|| =  0.0
In [58]:
plot_results(f_cormick, xi, yi, dim = [np.linspace(-45, 8, 1000), np.linspace(-45, 8, 1000)], method = '3')

Pelo Método BFGS, o mínimo da função foi encontrado com apenas 8 iterações, o que não ocorreu no Método de Newton, que resultou em um valor diferente do esperado. Se comparado com o Método do Gradiente, o número de iterações foi muito menor (46 iterações foram necessárias no primeiro Método), o que teoricamente faz sentido uma vez que o tipo de convergência é melhor se comparado com o Método do Gradiente.

Exemplo 4 - Função Colville¶

Testando agora na Função Colville:

In [59]:
f_colville_grad = nd.Gradient(f_colville)
In [60]:
chute = [3, 1, 5, 4]
xi, yi, grad_norm = BFGS(f_colville, f_colville_grad, chute = chute, tol=1e-8)  
df4_3 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'x3':xi[:,2], 'x4':xi[:,3], 'f(x1,x2,x3,x4)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'x3', 'x4', 'f(x1,x2,x3,x4)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Chute inicial: y = 46200.9, x = [3, 1, 5, 4] 


*************** Resultados ***************

Número de iterações:  77

***** Primeira iteração *****

f([3 1 5 4]) = 46200.9

||∇f([3 1 5 4])|| =  39215.92871169


***** Última iteração: 77° *****

f([1. 1. 1. 1.]) = 0.0

||∇f([1. 1. 1. 1.])|| =  0.0
In [61]:
df4_3
Out[61]:
x1 x2 x3 x4 f(x1,x2,x3,x4) ||∇f||
0 3.00000000 1.00000000 5.00000000 4.00000000 46200.90000000 39215.92871169
1 0.65527344 1.37612305 -4.23046875 4.90805664 15485.72247885 19923.11661312
2 -2.34055222 1.70232693 -3.63569501 7.61445932 4823.45883370 8221.27443435
3 -1.23940941 3.33747143 -3.66775263 11.54402139 2345.12224793 2751.91227485
4 -0.40102996 0.00289211 -3.86812964 13.23333383 1577.26670286 2426.65070686
... ... ... ... ... ... ...
72 0.99997852 0.99995828 1.00001805 1.00003609 0.00000000 0.00056475
73 0.99999924 0.99999858 1.00000052 1.00000093 0.00000000 0.00006695
74 1.00000000 1.00000001 0.99999999 0.99999998 0.00000000 0.00000135
75 1.00000000 1.00000000 1.00000000 1.00000000 0.00000000 0.00000001
76 1.00000000 1.00000000 1.00000000 1.00000000 0.00000000 0.00000000

77 rows × 6 columns

O mínimo da função foi encontrado com $78$ passos, mais que no Método de Newton.

Área de Testes¶

Nessa parte de testes, o primeiro passo é escrever a função que deseja minimizar. Para isso, siga:

$$ $$
f = lambda x: x[0]**2 + x[1]**2 + 2*x[0]*x[1] + x[2]
$$ $$

Ou

$$ $$
def f(x):
                                    return x[0]**2 + x[1]**2 + 2*x[0]*x[1] + x[2]
$$ $$

Para minimizar a função $f(x, y, z) = x^{2} + y^{2} + 2xy + z$, por exemplo.

$$ $$

Depois escolha o método a ser utilizado com um número inteiro (1, 2 ou 3):

  1. Método do Gradiente = 1

  2. Método de Newton = 2

  3. Método BFGS = 3

$$ $$

Finalmente, escolha o chute inicial (uma lista com a mesma dimensão da função), a constante de decréscimo de Armijo $\sigma$ (float $> 0$) e a tolerância (1e-8 para $10^{-8}$), por exemplo.

EXEMPLO¶

Função Objetivo:

$$ $$\begin{equation} f(\mathbf{x}) = (x_{1} + 2x_{2} - 7)^{2} + (2x_{1} + x_{2} - 5)^{2} + e^{-2y} \end{equation}$$ $$

Entrada da forma, para o Método do Gradiente:

$$ $$
f = lambda x: (x[0] + 2*x[1] - 7)**2 + (2*x[0] + x[1] - 5)**2 + np.exp(-2*x[1])
                            minimiza_f(f, metodo = 1, chute = [8, 15], sigma = 0.02, tol = 1e-10)
$$ $$

Cujo resultado está abaixo:

In [62]:
f = lambda x: (x[0] + 2*x[1] - 7)**2 + (2*x[0] + x[1] - 5)**2 + np.exp(-2*x[1])

metodo = 2

chute = [8, 15]

sigma = 0.5

tol = 1e-10

minimiza_f(f = f, metodo = metodo, chute = chute, sigma = sigma, tol = tol)
Método de Newton com condição de Armijo 



*************** Resultados ***************

Número de iterações:  24

***** Primeira iteração *****

f([ 8 15]) = 1637.0

||∇f([ 8 15])|| =  241.93387526


***** Última iteração: 24° *****

f([0.99890135 3.00137331]) = 0.00247535

||∇f([0.99890135 3.00137331])|| =  0.0





Referências¶

$$ $$

[1] SURJANOVIC, S., BINGHAM, D., Simon Fraser University. Virtual Library of Simulation Experiments: Test Functions and Datasets. [https://www.sfu.ca/~ssurjano/optimization.html]

$$ $$

[2] Test functions for optimization, Wikipedia. [https://en.wikipedia.org/wiki/Test_functions_for_optimization]

$$ $$

[3] Real Valued Test Functions, HeuristicLab. [https://dev.heuristiclab.com/trac.fcgi/wiki/Documentation/Reference/Test%20Functions] $$ $$

[4] Broyden–Fletcher–Goldfarb–Shanno algorithm, Wikipedia.[https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm]

$$ $$

[5] FRIEDLANDER, A. Elementos de Programação Não-Linear. Universidade Estadual de Campinas.[https://www.ime.unicamp.br/~friedlan/livro.pdf]